rm(list=ls())
library(tidyverse)
library(AER)
library(grid)
library(ggmap)
library(RColorBrewer)
library(ggpubr)
library(xtable)
library(stargazer)

##### PLOT LOCATION OF SKELETAL RECORDS IN EUROPE (FIG 1) ##### 
sites <- read.csv("SiteSkeletal_Euro.csv")

# Separate Iceland/Greenland from rest of Europe
sites_Europe <- filter(sites,
                       country != "Iceland" &
                         country != "Greenland")
sites_IceGreen <- filter(sites,
                         country %in% c("Iceland",
                                        "Greenland"))

# set min and max longitude/latitude for Europe map
min.lon.Eu <- min(sites_Europe$longitude, na.rm = T) - 1
max.lon.Eu <- max(sites_Europe$longitude, na.rm = T) + 1
min.lat.Eu <- min(sites_Europe$latitude, na.rm = T) - 1
max.lat.Eu <- max(sites_Europe$latitude, na.rm = T) + 1

# create base Europe map (w/o labels)
europe.p <- get_stamenmap(bbox = c(left = min.lon.Eu, 
                                   bottom = min.lat.Eu, 
                                   right = max.lon.Eu, 
                                   top = max.lat.Eu),
                          zoom = 4,
                          maptype = "toner-background") 
mapatt <- attributes(europe.p)
map_transparent <- matrix(adjustcolor(europe.p, alpha.f = 0.3), 
                          nrow = nrow(europe.p))
attributes(map_transparent) <- mapatt
europe.p <- ggmap(map_transparent)

bones.eu.p <- europe.p + 
  geom_point(data = sites_Europe, aes(x = longitude, y = latitude), 
             size = 2,
             shape = 10,
             alpha = 0.7,
             color = brewer.pal(3, "Dark2")[2]) + # add bones sites
  xlab("Longitude") +
  ylab("Latitude") +
  theme(panel.grid = element_blank(), 
        panel.background = element_rect(fill = 'white'))

# set min and max longitude/latitude for Iceland/Greenland map
min.lon.IG <- min(sites_IceGreen$longitude, na.rm = T) - 10
max.lon.IG <- max(sites_IceGreen$longitude, na.rm = T) + 8
min.lat.IG <- min(sites_IceGreen$latitude, na.rm = T) - 7
max.lat.IG <- max(sites_IceGreen$latitude, na.rm = T) + 10

# create base Iceland/Greenland map (w/o labels)
IG.p <- get_stamenmap(bbox = c(left = min.lon.IG, 
                               bottom = min.lat.IG, 
                               right = max.lon.IG, 
                               top = max.lat.IG),
                      zoom = 4,
                      maptype = "toner-background") 
mapatt <- attributes(IG.p)
map_transparent <- matrix(adjustcolor(IG.p, alpha.f = 0.3), 
                          nrow = nrow(IG.p))
attributes(map_transparent) <- mapatt
IG.p <- ggmap(map_transparent)

bones.IG.p <- IG.p + 
  geom_point(data = sites_IceGreen, aes(x = longitude, y = latitude), 
             size = 2,
             shape = 10,
             alpha = 0.7,
             color = brewer.pal(3, "Dark2")[2]) + # add bones sites
  xlab("Longitude") +
  ylab("Latitude") +
  theme(panel.grid = element_blank(), 
        panel.background = element_rect(fill = 'white'))

###### Export maps (Figure 1) ###### 
pdf("plots/bones_maps.pdf", width = 7, height = 3.5)
ggarrange(bones.eu.p, bones.IG.p)
dev.off()

##### PLOT REGIONAL DISTRIBUTION OF HISTORICAL PRO-FEMALE BIAS (FIG S1b.1) ##### 
#sites$settType3 <- ifelse(sites$settType2 %in% c("Urban, large", "Urban, small"), "Urban",
#                          "Rural")
region.LEH <- sites %>% group_by(region, settType2) %>%
  summarize(nMales = sum(nMales, na.rm = T),
            nFemales = sum(nFemales, na.rm = T),
            nMalesAffected = sum(nMalesAffected, na.rm = T),
            nFemalesAffected = sum(nFemalesAffected, na.rm = T)) %>%
  mutate(lehDiff = nMalesAffected/nMales - nFemalesAffected/nFemales)
region.LEH <- filter(region.LEH, settType2 != "")

region.LEH.p <- ggplot(region.LEH, aes(x = region, y = lehDiff, fill = settType2)) +
  geom_bar(stat="identity", position=position_dodge(),
           color="black") +
  scale_x_discrete(name = "Region",
                   labels = c("Central and\nSoutheastern\nEurope",
                              "Mediterranean",
                              "Northeastern\nEurope",
                              "Northern\nEurope",
                              "Northwestern\nEurope")) +
  scale_fill_brewer(name = "Settlement\nType",
                    labels = c("Rural", "Urban (large)",
                               "Urban (small)"),
                    palette=6) +
  scale_y_continuous(name = "Historical Pro-Female Bias") +
  theme_bw() +
  theme(text=element_text(size=6))

###### Export regional distribution plot (Figure S1b.1) ###### 
pdf("plots/region.LEH.p.pdf", width = 4, height = 2)
region.LEH.p
dev.off()

##### DESCRIPTIVE STATISTICS (TABLES S1f.1 & S1f.2) ##### 
###### Load EVS data ######
rm(list=ls())
load("DataAnalysis_Euro.RData")

###### Descriptive Statistics for EVS Wave 4 (Table S1f.1) ######
table(wave4.Church$settType2)
wave4.Church$Rural <- ifelse(wave4.Church$settType2 == "Rural", 1,
                          ifelse(wave4.Church$settType2 %in% c("Urban, large", "Urban, small"), 0, NA))
table(wave4.Church$Rural)
wave4.Church$Urban_small <- ifelse(wave4.Church$settType2 == "Urban, small", 1,
                                ifelse(wave4.Church$settType2 %in% c("Urban, large", "Rural"), 0, NA))
table(wave4.Church$Urban_small)
wave4.Church$Urban_large <- ifelse(wave4.Church$settType2 == "Urban, large", 1,
                                ifelse(wave4.Church$settType2 %in% c("Urban, small", "Rural"), 0, NA))
table(wave4.Church$specialty)
wave4.Church$CraftArtisan <- ifelse(wave4.Church$specialty == "Craft or artisan workers", 1,0)
table(wave4.Church$CraftArtisan)
wave4.Church$Farming <- ifelse(wave4.Church$specialty == "Farming community", 1,0)
table(wave4.Church$Farming)
wave4.Church$Hospital <- ifelse(wave4.Church$specialty == "Hospital population", 1,0)
table(wave4.Church$Hospital)
wave4.Church$Other <- ifelse(wave4.Church$specialty == "Other", 1,0)
table(wave4.Church$Other)
wave4.Church$Religious <- ifelse(wave4.Church$specialty == "Religious order", 1,0)
table(wave4.Church$Religious)

sum_stats_w4 <- wave4.Church %>%
  dplyr::select(allDVs,stereotype,malePreference,
                womenFulfill.x,workMom.x,
                housewife.x,womenWant.x,
                fatherSuited.x,workMomSuffer.x,
                spouses.x,womenJob.x,menWork.x,
                lehDiff,plagueSeverity,
                lowElevation,topography.x,
                postChristian,medieval,postMedieval,
                Rural,Urban_small,Urban_large,
                suitability,CraftArtisan,Farming,
                Hospital,Religious,Other,
                sex,age,unemployed,incomeDecile,
                educ,married,religiosity,
                ChurchExp) %>%
  psych::describe() %>%
  data.frame() # summary statistics
sum_stats_w4 <- sum_stats_w4[,c(2:5, 8:9,11)]
rownames(sum_stats_w4) <- c("Modern Pro-Female Bias",
                            "Gender Roles",
                            "Pro-Female Preference",
                            "Women Child",
                            "Working Mom",
                            "Housewife",
                            "Women Want",
                            "Father Child",
                            "Child Suffer",
                            "Dual Income",
                            "Independent Women",
                            "Men Work",
                            "Historical Pro-Female Bias",
                            "Black Death Severity",
                            "Low Elevation",
                            "Topography",
                            "Post Christian",
                            "Medieval",
                            "Post-Medieval",
                            "Rural",
                            "Urban, small",
                            "Urban, large",
                            "Crop Suitability",
                            "Specialty: Craft Artisan",
                            "Specialty: Farming",
                            "Specialty: Hospital",
                            "Specialty: Religious Order",
                            "Specialty: Other",
                            "Sex",
                            "Age",
                            "Unemployed",
                            "Income",
                            "Education",
                            "Married",
                            "Religiosity",
                            "Western Church Exposure")
xtable(sum_stats_w4)

###### Descriptive Statistics for EVS Wave 5 (Table S1f.2) ######
table(wave5.Church$settType2)
wave5.Church$Rural <- ifelse(wave5.Church$settType2 == "Rural", 1,
                          ifelse(wave5.Church$settType2 %in% c("Urban, large", "Urban, small"), 0, NA))
table(wave5.Church$Rural)
wave5.Church$Urban_small <- ifelse(wave5.Church$settType2 == "Urban, small", 1,
                                ifelse(wave5.Church$settType2 %in% c("Urban, large", "Rural"), 0, NA))
table(wave5.Church$Urban_small)
wave5.Church$Urban_large <- ifelse(wave5.Church$settType2 == "Urban, large", 1,
                                ifelse(wave5.Church$settType2 %in% c("Urban, small", "Rural"), 0, NA))
table(wave5.Church$specialty)
wave5.Church$CraftArtisan <- ifelse(wave5.Church$specialty == "Craft or artisan workers", 1,0)
table(wave5.Church$CraftArtisan)
wave5.Church$Farming <- ifelse(wave5.Church$specialty == "Farming community", 1,0)
table(wave5.Church$Farming)
#wave5.new$Hospital <- ifelse(wave5.new$specialty == "Hospital population", 1,0)
#table(wave5.new$Hospital) # no obs living in previous hospital site
wave5.Church$Other <- ifelse(wave5.Church$specialty == "Other", 1,0)
table(wave5.Church$Other)
wave5.Church$Religious <- ifelse(wave5.Church$specialty == "Religious order", 1,0)
table(wave5.Church$Religious)

sum_stats_w5 <- wave5.Church %>%
  dplyr::select(allDVs,stereotype,malePreference,
                womenWant.x,workMomSuffer.x,
                menWork.x,menLead.x,
                menEduc.x,menExec.x,
                womenRights.x,
                lehDiff,plagueSeverity,
                lowElevation,topography.x,
                postChristian,medieval,postMedieval,
                Rural,Urban_small,Urban_large,
                suitability,CraftArtisan,Farming,
                Religious,Other,
                sex,age,unemployed,incomeDecile,
                educ,married,religiosity,
                ChurchExp) %>%
  psych::describe() %>%
  data.frame() # summary statistics
sum_stats_w5 <- sum_stats_w5[,c(2:5, 8:9,11)]
rownames(sum_stats_w5) <- c("Modern Pro-Female Bias",
                            "Gender Roles",
                            "Pro-Female Preference",
                            "Women Want",
                            "Child Suffer",
                            "Men Work",
                            "Men Leader",
                            "Men Education",
                            "Men Executives",
                            "Women Rights",
                            "Historical Pro-Female Bias",
                            "Black Death Severity",
                            "Low Elevation",
                            "Topography",
                            "Post Christian",
                            "Medieval",
                            "Post-Medieval",
                            "Rural",
                            "Urban, small",
                            "Urban, large",
                            "Crop Suitability",
                            "Specialty: Craft Artisan",
                            "Specialty: Farming",
                            "Specialty: Religious Order",
                            "Specialty: Other",
                            "Sex",
                            "Age",
                            "Unemployed",
                            "Income",
                            "Education",
                            "Married",
                            "Religiosity",
                            "Western Church Exposure")
xtable(sum_stats_w5)

###### Within and between country variances (Table S1f.3) ######
with.bet.lehDiff.w4 <- lm(lehDiff ~ factor(country),
                          data = wave4.new)
anova(with.bet.lehDiff.w4)["Residuals", "Mean Sq"] # within-country
anova(with.bet.lehDiff.w4)["factor(country)", "Mean Sq"] # between-country
anova(with.bet.lehDiff.w4)["Residuals", "Mean Sq"]/anova(with.bet.lehDiff.w4)["factor(country)", "Mean Sq"]

with.bet.lehDiff.w5 <- lm(lehDiff ~ factor(country),
                          data = wave5.new)
anova(with.bet.lehDiff.w5)["Residuals", "Mean Sq"] # within-country
anova(with.bet.lehDiff.w5)["factor(country)", "Mean Sq"] # between-country
anova(with.bet.lehDiff.w5)["Residuals", "Mean Sq"]/anova(with.bet.lehDiff.w5)["factor(country)", "Mean Sq"]

with.bet.allDVs.w4 <- lm(allDVs ~ factor(country),
                          data = wave4.new)
anova(with.bet.allDVs.w4)["Residuals", "Mean Sq"] # within-country
anova(with.bet.allDVs.w4)["factor(country)", "Mean Sq"] # between-country
anova(with.bet.allDVs.w4)["Residuals", "Mean Sq"]/anova(with.bet.allDVs.w4)["factor(country)", "Mean Sq"]

with.bet.allDVs.w5 <- lm(allDVs ~ factor(country),
                         data = wave5.new)
anova(with.bet.allDVs.w5)["Residuals", "Mean Sq"] # within-country
anova(with.bet.allDVs.w5)["factor(country)", "Mean Sq"] # between-country
anova(with.bet.allDVs.w5)["Residuals", "Mean Sq"]/anova(with.bet.allDVs.w5)["factor(country)", "Mean Sq"]

##### BASELINE EVS RESULTS USING SEQUENTIAL-G ##### 
# create function to obtain bootstrapped SEs
boot.SE <- function(data, # specify data to be analyzed on
                    boots, # number of iterations 
                    x, # specify predictor variable (character variable)
                    y, # specify outcome variable (character variable)
                    FE = "+ as.factor(country)", # specify FE (character)
                    pretreat.formula, # specify pretreat formula (RHS only, in character class)
                    posttreat.formula) # specify posttreat formula (RHS only, in character class)
{ 
  posttreat <- unlist(str_split(posttreat.formula, "\\+")) # create vector of posttreat variable names
  fs_main <- lm(as.formula(paste(y, 
                                 "~",
                                 x,
                                 FE,
                                 "+",
                                 pretreat.formula,
                                 "+",
                                 posttreat.formula)),
                data = data,
                weights=1/data$nObs, 
                na.action=na.exclude)
  demed.y_main <- data[,y] - as.matrix(data[,posttreat]) %*% coef(fs_main)[posttreat] # get demediated outcome
  data$demed.y_main <- as.numeric(unlist(demed.y_main)) # add demed y to databoot
  ss_main <- lm(as.formula(paste("demed.y_main",
                                 "~",
                                 x,
                                 "+",
                                 pretreat.formula)), # run second stage regression
                data = data,
                weights=1/data$nObs, 
                na.action=na.exclude)
  gboot <- matrix(NA, nrow=boots, ncol= length(names(coef(ss_main)))) # empty matrix to store bootstrapped values
  set.seed(12345)
  for (b in 1:boots){
    databoot <- data[sample(1:nrow(data), replace = TRUE),]
    fs <- lm(as.formula(paste(y, 
                              "~",
                              x,
                              FE,
                              "+",
                              pretreat.formula,
                              "+",
                              posttreat.formula)), # run first stage regression
             data = databoot,
             weights=1/databoot$nObs, 
             na.action=na.exclude)
    demed.y <- databoot[,y] - as.matrix(databoot[,posttreat]) %*% coef(fs)[posttreat] # get demediated outcome
    databoot$demed.y <- as.numeric(unlist(demed.y)) # add demed y to databoot
    ss <- lm(as.formula(paste("demed.y",
                              "~",
                              x,
                              "+",
                              pretreat.formula)), # run second stage regression
             data = databoot,
             weights=1/databoot$nObs, 
             na.action=na.exclude)
    gboot[b,] <- coef(ss) # store bootstrapped coefs
  }
  colnames(gboot) <- names(coef(ss))
  return(list(model = ss_main,
              se = apply(gboot,2,sd))) # return bootstrapped SEs
}

# set baseline model specifications (make sure no spaces between variables)
pretreat.model <- "lowElevation+as.factor(topography.x)+postChristian+medieval+postMedieval+as.factor(settType2)"
posttreat.model <- "sex+age+unemployed+incomeDecile+married+religiosity+educ"

###### Baseline models with sequential-g estimation ###### 
m.w4.allDVs <- boot.SE(data = wave4.new,
                       boots = 1000,
                       x = "lehDiff",
                       y = "allDVs",
                       pretreat.formula = pretreat.model,
                       posttreat.formula = posttreat.model)

m.w5.allDVs <- boot.SE(data = wave5.new,
                       boots = 1000,
                       x = "lehDiff",
                       y = "allDVs",
                       pretreat.formula = pretreat.model,
                       posttreat.formula = posttreat.model)

m.w4.malePref <- boot.SE(data = wave4.new,
                         boots = 1000,
                         x = "lehDiff",
                         y = "malePreference",
                         pretreat.formula = pretreat.model,
                         posttreat.formula = posttreat.model)

m.w5.malePref <- boot.SE(data = wave5.new,
                         boots = 1000,
                         x = "lehDiff",
                         y = "malePreference",
                         pretreat.formula = pretreat.model,
                         posttreat.formula = posttreat.model)

m.w4.stereotype <- boot.SE(data = wave4.new,
                           boots = 1000,
                           x = "lehDiff",
                           y = "stereotype",
                           pretreat.formula = pretreat.model,
                           posttreat.formula = posttreat.model)

m.w5.stereotype <- boot.SE(data = wave5.new,
                           boots = 1000,
                           x = "lehDiff",
                           y = "stereotype",
                           pretreat.formula = pretreat.model,
                           posttreat.formula = posttreat.model)

# Export regression tables for SI (Table S2b.1)
stargazer(m.w4.allDVs$model,
          m.w4.stereotype$model,
          m.w4.malePref$model,
          style = "ajps",
          se = list(m.w4.allDVs$se,
                    m.w4.stereotype$se,
                    m.w4.malePref$se),
          title = "Wave 4 Full Results",
          model.numbers = FALSE,
          digits = 3,
          omit.stat = c("f", "ser", "aic", "wald", 
                        "rsq", "adj.rsq"),
          star.char = c("*", "**"),
          star.cutoffs = c(0.05, 0.01))

stargazer(m.w5.allDVs$model,
          m.w5.stereotype$model,
          m.w5.malePref$model,
          style = "ajps",
          se = list(m.w5.allDVs$se,
                    m.w5.stereotype$se,
                    m.w5.malePref$se),
          title = "Wave 5 Full Results",
          model.numbers = FALSE,
          digits = 3,
          omit.stat = c("f", "ser", "aic", "wald", 
                        "rsq", "adj.rsq"),
          star.char = c("*", "**"),
          star.cutoffs = c(0.05, 0.01))

###### Plot coefficient estimates for main paper (Figure 2) ###### 
coef_out.df <- data.frame(waves = rep(c("Wave 4", "Wave 5"),
                                      each = 3),
                          outcome = rep(c("allDVs",
                                          "malePreference",
                                          "stereotype"),
                                        2),
                          coef = NA,
                          se = NA,
                          lb = NA,
                          ub = NA)

coef_out.df$outcome <- factor(coef_out.df$outcome,
                              levels = rev(c("allDVs",
                                             "stereotype",
                                             "malePreference")))

coef_out.df$coef[coef_out.df$waves == "Wave 4" & 
                   coef_out.df$outcome == "allDVs"] <- coef(m.w4.allDVs$model)["lehDiff"]
coef_out.df$coef[coef_out.df$waves == "Wave 5" & 
                   coef_out.df$outcome == "allDVs"] <- coef(m.w5.allDVs$model)["lehDiff"]
coef_out.df$coef[coef_out.df$waves == "Wave 4" & 
                   coef_out.df$outcome == "stereotype"] <- coef(m.w4.stereotype$model)["lehDiff"]
coef_out.df$coef[coef_out.df$waves == "Wave 5" & 
                   coef_out.df$outcome == "stereotype"] <- coef(m.w5.stereotype$model)["lehDiff"]
coef_out.df$coef[coef_out.df$waves == "Wave 4" & 
                   coef_out.df$outcome == "malePreference"] <- coef(m.w4.malePref$model)["lehDiff"]
coef_out.df$coef[coef_out.df$waves == "Wave 5" & 
                   coef_out.df$outcome == "malePreference"] <- coef(m.w5.malePref$model)["lehDiff"]

coef_out.df$se[coef_out.df$waves == "Wave 4" & 
                 coef_out.df$outcome == "allDVs"] <- m.w4.allDVs$se["lehDiff"]
coef_out.df$se[coef_out.df$waves == "Wave 5" & 
                 coef_out.df$outcome == "allDVs"] <- m.w5.allDVs$se["lehDiff"]
coef_out.df$se[coef_out.df$waves == "Wave 4" & 
                 coef_out.df$outcome == "stereotype"] <- m.w4.stereotype$se["lehDiff"]
coef_out.df$se[coef_out.df$waves == "Wave 5" & 
                 coef_out.df$outcome == "stereotype"] <- m.w5.stereotype$se["lehDiff"]
coef_out.df$se[coef_out.df$waves == "Wave 4" & 
                 coef_out.df$outcome == "malePreference"] <- m.w4.malePref$se["lehDiff"]
coef_out.df$se[coef_out.df$waves == "Wave 5" & 
                 coef_out.df$outcome == "malePreference"] <- m.w5.malePref$se["lehDiff"]

coef_out.df$lb <- coef_out.df$coef - 1.96*coef_out.df$se
coef_out.df$ub <- coef_out.df$coef + 1.96*coef_out.df$se

text <- textGrob("Sub-indices", rot = 90,
                 gp=gpar(fontsize=5))

coef.out.p <- ggplot(coef_out.df, aes(x = coef, y = outcome)) +
  geom_point(size = 0.8) +
  geom_errorbarh(aes(xmin = lb, xmax = ub),
                 height = 0) +
  facet_wrap("waves", scale="fixed") +
  geom_vline(xintercept = 0,
             linetype = "dashed") +
  scale_y_discrete(name = "",
                   labels = c("Pro-Female\nPreference",
                              "Gender Roles",
                              "Modern\nPro-Female Bias")) +
  scale_x_continuous(name = "Effect of Historical Pro-Female Bias",
                     limits = c(-0.1,0.5)) +
  coord_cartesian(ylim = c(1, 3), clip="off") +
  annotation_custom(text, xmin = -1.3, xmax = -1.3, ymin = 1.5, ymax = 1.5) +
  annotation_custom(linesGrob(), xmin = -1.25, xmax = -1.25, ymin = 2, ymax = 1) +
  annotation_custom(linesGrob(), xmin = -1.25, xmax = -1.23, ymin = 2, ymax = 2) +
  annotation_custom(linesGrob(), xmin = -1.25, xmax = -1.23, ymin = 1, ymax = 1) +
  theme_bw() +
  theme(text=element_text(size=6),
        plot.margin = unit(c(1,1,1,1), "lines"))

# Export coef plot for main EVS analysis (Figure 2)
pdf("plots/mainEVS.pdf", width = 3, height = 2)
coef.out.p
dev.off()

##### ROBUSTNESS CHECKS 1: OLS WITH PRETREATMENT COVARIATES ONLY ##### 
###### Wave 4: DV = allDVs ###### 
m1.allDVs.w4 <- lm(allDVs ~ lehDiff+
                     lowElevation+as.factor(topography.x)+postChristian+medieval+postMedieval+as.factor(settType2),
                   data = wave4.new,
                   weights = wave4.new$nObs)
se.m1.allDVs.w4 <- sqrt(diag(vcovHC(m1.allDVs.w4, type = "HC0")))
coeftest(m1.allDVs.w4, vcovHC(m1.allDVs.w4, type = "HC0"))

###### Wave 4: DV = malePreference ###### 
m1.malePreference.w4 <- lm(malePreference ~ lehDiff+
                             lowElevation+as.factor(topography.x)+postChristian+medieval+postMedieval+as.factor(settType2),
                           data = wave4.new,
                           weights = wave4.new$nObs)
se.m1.malePreference.w4 <- sqrt(diag(vcovHC(m1.malePreference.w4, type = "HC0")))
coeftest(m1.malePreference.w4, vcovHC(m1.malePreference.w4,type = "HC0"))

###### Wave 4: DV = stereotype ###### 
m1.stereotype.w4 <- lm(stereotype ~ lehDiff+
                         lowElevation+as.factor(topography.x)+postChristian+medieval+postMedieval+as.factor(settType2),
                       data = wave4.new,
                       weights = wave4.new$nObs)
se.m1.stereotype.w4 <- sqrt(diag(vcovHC(m1.stereotype.w4, type = "HC0")))
coeftest(m1.stereotype.w4, vcovHC(m1.stereotype.w4,
                                  type = "HC0"))

###### Wave 5: DV = allDVs ###### 
m1.allDVs.w5 <- lm(allDVs ~ lehDiff+
                     lowElevation+as.factor(topography.x)+postChristian+medieval+postMedieval+as.factor(settType2),
                   data = wave5.new,
                   weights = wave5.new$nObs)
se.m1.allDVs.w5 <- sqrt(diag(vcovHC(m1.allDVs.w5, type = "HC0")))
coeftest(m1.allDVs.w5, vcovHC(m1.allDVs.w5,
                              type = "HC0"))

###### Wave 5: DV = malePreference ###### 
m1.malePreference.w5 <- lm(malePreference ~ lehDiff+
                             lowElevation+as.factor(topography.x)+postChristian+medieval+postMedieval+as.factor(settType2),
                           data = wave5.new,
                           weights = wave5.new$nObs)
se.m1.malePreference.w5 <- sqrt(diag(vcovHC(m1.malePreference.w5, type = "HC0")))
coeftest(m1.malePreference.w5, vcovHC(m1.malePreference.w5,
                                      type = "HC0"))

###### Wave 5: DV = stereotype ###### 
m1.stereotype.w5 <- lm(stereotype ~ lehDiff+
                         lowElevation+as.factor(topography.x)+postChristian+medieval+postMedieval+as.factor(settType2),
                       data = wave5.new,
                       weights = wave5.new$nObs)
se.m1.stereotype.w5 <- sqrt(diag(vcovHC(m1.stereotype.w5, type = "HC0")))
coeftest(m1.stereotype.w5, vcovHC(m1.stereotype.w5,
                                  type = "HC0"))

# Export regression tables for SI (Table S2d.1)
stargazer(m1.allDVs.w4,
          m1.stereotype.w4,
          m1.malePreference.w4,
          style = "ajps",
          se = list(se.m1.allDVs.w4,
                    se.m1.stereotype.w4,
                    se.m1.malePreference.w4),
          title = "Wave 4 Full Results (OLS Regression)",
          model.numbers = FALSE,
          digits = 3,
          omit.stat = c("f", "ser", "aic", "wald", 
                        "rsq", "adj.rsq"),
          star.char = c("*", "**"),
          star.cutoffs = c(0.05, 0.01))

stargazer(m1.allDVs.w5,
          m1.stereotype.w5,
          m1.malePreference.w5,
          style = "ajps",
          se = list(se.m1.allDVs.w5,
                    se.m1.stereotype.w5,
                    se.m1.malePreference.w5),
          title = "Wave 5 Full Results (OLS Regression)",
          model.numbers = FALSE,
          digits = 3,
          omit.stat = c("f", "ser", "aic", "wald", 
                        "rsq", "adj.rsq"),
          star.char = c("*", "**"),
          star.cutoffs = c(0.05, 0.01))

##### ROBUSTNESS CHECKS 2: EVS-ITEM REGRESSION USING SEQUENTIAL-G ##### 
coef_out_items.df <- data.frame(waves = rep(c("wave4.new",
                                              "wave5.new"),
                                            each = 13),
                                outcome = rep(c("menWork.x","womenFulfill.x",
                                                "housewife.x","spouses.x",
                                                "menLead.x","menEduc.x",
                                                "workMom.x","workMomSuffer.x",
                                                "womenWant.x","womenJob.x",
                                                "fatherSuited.x","menExec.x",
                                                "womenRights.x"),
                                              2),
                                coef = NA,
                                lb = NA,
                                ub = NA)

for (i in 1:nrow(coef_out_items.df)){
  data <- eval(parse(text = coef_out_items.df$waves[i]))
  if (sum(!is.na(data[,coef_out_items.df$outcome[i]])) > 2000){
    model <- boot.SE(data = data,
                     boots = 1000,
                     x = "lehDiff",
                     y = coef_out_items.df$outcome[i],
                     pretreat.formula = pretreat.model,
                     posttreat.formula = posttreat.model)
    coef_out_items.df$coef[i] <- coef(model$model)["lehDiff"]
    coef_out_items.df$lb[i] <- coef(model$model)["lehDiff"] - model$se["lehDiff"] * 1.96
    coef_out_items.df$ub[i] <- coef(model$model)["lehDiff"] + model$se["lehDiff"] * 1.96
  } else {
    coef_out_items.df$coef[i] <- NA
    coef_out_items.df$lb[i] <- NA
    coef_out_items.df$ub[i] <- NA
  }
}

coef_out_items.df$labels <- rep(c("Men Work",
                                  "Women Child",
                                  "Housewife",
                                  "Dual Income",
                                  "Men Leader",
                                  "Men Education",
                                  "Working Mom",
                                  "Child Suffer",
                                  "Women Want",
                                  "Independent\nWomen",
                                  "Father Child",
                                  "Men Executives",
                                  "Women Rights"),2)
coef_out_items.df$group <- rep(c("Pro-Female\nPreference",
                                 rep("Gender\nRoles",3),
                                 rep("Pro-Female\nPreference",2),
                                 rep("Gender\nRoles",5),
                                 "Pro-Female\nTreatment",
                                 "Others"),2)
coef_out_items.df$group <- factor(coef_out_items.df$group,
                                  levels = c("Gender\nRoles",
                                             "Pro-Female\nPreference",
                                             "Others"))

coef_out_items.df.clean <- na.omit(coef_out_items.df)
coef_out_items.p <- ggplot(coef_out_items.df.clean, 
                           aes(x = coef, y = labels, 
                               shape = group,
                               linetype = group)) +
  geom_point() +
  geom_errorbarh(aes(xmin = lb, xmax = ub),
                 height = 0) + 
  facet_wrap("waves",
             scales = "free_y",
             labeller = as_labeller(c("wave4.new" = "Wave 4",
                                      "wave5.new" = "Wave 5"))) +
  geom_vline(xintercept = 0,
             linetype = "dashed") +
  scale_y_discrete(name = "") +
  scale_x_continuous(name = "Effect of Historical Pro-Female Bias") +
  scale_shape_discrete(name = "Sub-indices") +
  scale_linetype_discrete(name = "Sub-indices") +
  theme_bw() +
  theme(text = element_text(size=6),
        legend.position="top")

# Export coef plot for individual gender item (Figure S2d.1)
pdf("plots/coef_out_items.pdf", width = 4, height = 3)
coef_out_items.p
dev.off()

##### ROBUSTNESS CHECKS 3: REMOVE RECENT SITES ##### 
###### Sequential-g models for DVs ###### 
m.w4.allDVs.recent <- boot.SE(data = filter(wave4.new, distanceFromNow >=200),
                              boots = 1000,
                              x = "lehDiff",
                              y = "allDVs",
                              pretreat.formula = pretreat.model,
                              posttreat.formula = posttreat.model)

m.w5.allDVs.recent <- boot.SE(data = filter(wave5.new, distanceFromNow >=200),
                              boots = 1000,
                              x = "lehDiff",
                              y = "allDVs",
                              pretreat.formula = pretreat.model,
                              posttreat.formula = posttreat.model)

m.w4.malePref.recent <- boot.SE(data = filter(wave4.new, distanceFromNow >=200),
                                boots = 1000,
                                x = "lehDiff",
                                y = "malePreference",
                                pretreat.formula = pretreat.model,
                                posttreat.formula = posttreat.model)

m.w5.malePref.recent <- boot.SE(data = filter(wave5.new, distanceFromNow >=200),
                                boots = 1000,
                                x = "lehDiff",
                                y = "malePreference",
                                pretreat.formula = pretreat.model,
                                posttreat.formula = posttreat.model)

m.w4.stereotype.recent <- boot.SE(data = filter(wave4.new, distanceFromNow >=200),
                                  boots = 1000,
                                  x = "lehDiff",
                                  y = "stereotype",
                                  pretreat.formula = pretreat.model,
                                  posttreat.formula = posttreat.model)

m.w5.stereotype.recent <- boot.SE(data = filter(wave5.new, distanceFromNow >=200),
                                  boots = 1000,
                                  x = "lehDiff",
                                  y = "stereotype",
                                  pretreat.formula = pretreat.model,
                                  posttreat.formula = posttreat.model)

# Export regression tables for SI (Table S2d.2)
stargazer(m.w4.allDVs.recent$model,
          m.w4.stereotype.recent$model,
          m.w4.malePref.recent$model,
          style = "ajps",
          se = list(m.w4.allDVs.recent$se,
                    m.w4.stereotype.recent$se,
                    m.w4.malePref.recent$se),
          title = "Wave 4 Results (Exclude Recent Sites)",
          model.numbers = FALSE,
          digits = 3,
          omit.stat = c("f", "ser", "aic", "wald", 
                        "rsq", "adj.rsq"),
          star.char = c("*", "**"),
          star.cutoffs = c(0.05, 0.01))

stargazer(m.w5.allDVs.recent$model,
          m.w5.stereotype.recent$model,
          m.w5.malePref.recent$model,
          style = "ajps",
          se = list(m.w5.allDVs.recent$se,
                    m.w5.stereotype.recent$se,
                    m.w5.malePref.recent$se),
          title = "Wave 4 and 5 Results (Exclude Recent Sites)",
          model.numbers = FALSE,
          digits = 3,
          omit.stat = c("f", "ser", "aic", "wald", 
                        "rsq", "adj.rsq"),
          star.char = c("*", "**"),
          star.cutoffs = c(0.05, 0.01))

##### ROBUSTNESS CHECKS 4: USE REGION FIXED EFFECTS ##### 
###### Sequential-g models for DVs (w region fixed effects) ###### 
m.w4.allDVs.coord <- boot.SE(data = wave4.new,
                             boots = 1000,
                             x = "lehDiff",
                             y = "allDVs",
                             FE = "+ as.factor(region)",
                             pretreat.formula = pretreat.model,
                             posttreat.formula = posttreat.model)

m.w5.allDVs.coord <- boot.SE(data = wave5.new,
                             boots = 1000,
                             x = "lehDiff",
                             y = "allDVs",
                             FE = "+ as.factor(region)",
                             pretreat.formula = pretreat.model,
                             posttreat.formula = posttreat.model)

m.w4.malePref.coord <- boot.SE(data = wave4.new,
                               boots = 1000,
                               x = "lehDiff",
                               y = "malePreference",
                               FE = "+ as.factor(region)",
                               pretreat.formula = pretreat.model,
                               posttreat.formula = posttreat.model)

m.w5.malePref.coord <- boot.SE(data = wave5.new,
                               boots = 1000,
                               x = "lehDiff",
                               y = "malePreference",
                               FE = "+ as.factor(region)",
                               pretreat.formula = pretreat.model,
                               posttreat.formula = posttreat.model)

m.w4.stereotype.coord <- boot.SE(data = wave4.new,
                                 boots = 1000,
                                 x = "lehDiff",
                                 y = "stereotype",
                                 FE = "+ as.factor(region)",
                                 pretreat.formula = pretreat.model,
                                 posttreat.formula = posttreat.model)

m.w5.stereotype.coord <- boot.SE(data = wave5.new,
                                 boots = 1000,
                                 x = "lehDiff",
                                 y = "stereotype",
                                 FE = "+ as.factor(region)",
                                 pretreat.formula = pretreat.model,
                                 posttreat.formula = posttreat.model)

# Export regression tables for SI (Table S2d.3)
stargazer(m.w4.allDVs.coord$model,
          m.w4.stereotype.coord$model,
          m.w4.malePref.coord$model,
          m.w5.allDVs.coord$model,
          m.w5.stereotype.coord$model,
          m.w5.malePref.coord$model,
          style = "ajps",
          se = list(m.w4.allDVs.coord$se,
                    m.w4.stereotype.coord$se,
                    m.w4.malePref.coord$se,
                    m.w5.allDVs.coord$se,
                    m.w5.stereotype.coord$se,
                    m.w5.malePref.coord$se),
          title = "Wave 4 and 5 Results (Use Region Fixed Effects)",
          model.numbers = FALSE,
          digits = 3,
          omit.stat = c("f", "ser", "aic", "wald", 
                        "rsq", "adj.rsq"),
          star.char = c("*", "**"),
          star.cutoffs = c(0.05, 0.01))

##### ROBUSTNESS CHECKS 5: ALTERNATIVE MEASURE OF CROP SUITABILITY ##### 
###### Set new pretreatment model (w alternative crop measure) ###### 
pretreat.model.crop <- "lowElevation+as.factor(topography.x)+suitability+postChristian+medieval+postMedieval+as.factor(settType2)"

###### Sequential-g models for DVs (w alternative crop measure) ###### 
m.w4.allDVs.crop <- boot.SE(data = wave4.new,
                            boots = 1000,
                            x = "lehDiff",
                            y = "allDVs",
                            pretreat.formula = pretreat.model.crop,
                            posttreat.formula = posttreat.model)

m.w5.allDVs.crop <- boot.SE(data = wave5.new,
                              boots = 1000,
                              x = "lehDiff",
                              y = "allDVs",
                              pretreat.formula = pretreat.model.crop,
                              posttreat.formula = posttreat.model)

m.w4.malePref.crop <- boot.SE(data = wave4.new,
                              boots = 1000,
                              x = "lehDiff",
                              y = "malePreference",
                              pretreat.formula = pretreat.model.crop,
                              posttreat.formula = posttreat.model)

m.w5.malePref.crop <- boot.SE(data = wave5.new,
                              boots = 1000,
                              x = "lehDiff",
                              y = "malePreference",
                              pretreat.formula = pretreat.model.crop,
                              posttreat.formula = posttreat.model)

m.w4.stereotype.crop <- boot.SE(data = wave4.new,
                                boots = 1000,
                                x = "lehDiff",
                                y = "stereotype",
                                pretreat.formula = pretreat.model.crop,
                                posttreat.formula = posttreat.model)

m.w5.stereotype.crop <- boot.SE(data = wave5.new,
                                boots = 1000,
                                x = "lehDiff",
                                y = "stereotype",
                                pretreat.formula = pretreat.model.crop,
                                posttreat.formula = posttreat.model)

# Export regression tables for SI (Table S2d.4)
stargazer(m.w4.allDVs.crop$model,
          m.w4.stereotype.crop$model,
          m.w4.malePref.crop$model,
          m.w5.allDVs.crop$model,
          m.w5.stereotype.crop$model,
          m.w5.malePref.crop$model,
          style = "ajps",
          se = list(m.w4.allDVs.crop$se,
                    m.w4.stereotype.crop$se,
                    m.w4.malePref.crop$se,
                    m.w5.allDVs.crop$se,
                    m.w5.stereotype.crop$se,
                    m.w5.malePref.crop$se),
          title = "Wave 4 and 5 Results (Alternative Measure of Crop Suitability)",
          model.numbers = FALSE,
          digits = 3,
          omit.stat = c("f", "ser", "aic", "wald", 
                        "rsq", "adj.rsq"),
          star.char = c("*", "**"),
          star.cutoffs = c(0.05, 0.01))

##### ROBUSTNESS CHECKS 6: ALTERNATIVE MEASURE OF SETTLEMENT TYPE ##### 
###### Set new pretreatment model (w alternative measure of settlement type) ###### 
pretreat.model.spec <- "lowElevation+as.factor(topography.x)+postChristian+medieval+postMedieval+as.factor(specialty)"

###### Sequential-g models for DVs (w alternative measure of settlement type) ###### 
m.w4.allDVs.spec <- boot.SE(data = wave4.new,
                            boots = 1000,
                            x = "lehDiff",
                            y = "allDVs",
                            pretreat.formula = pretreat.model.spec,
                            posttreat.formula = posttreat.model)

m.w5.allDVs.spec <- boot.SE(data = wave5.new,
                            boots = 1000,
                            x = "lehDiff",
                            y = "allDVs",
                            pretreat.formula = pretreat.model.spec,
                            posttreat.formula = posttreat.model)

m.w4.malePref.spec <- boot.SE(data = wave4.new,
                              boots = 1000,
                              x = "lehDiff",
                              y = "malePreference",
                              pretreat.formula = pretreat.model.spec,
                              posttreat.formula = posttreat.model)

m.w5.malePref.spec <- boot.SE(data = wave5.new,
                              boots = 1000,
                              x = "lehDiff",
                              y = "malePreference",
                              pretreat.formula = pretreat.model.spec,
                              posttreat.formula = posttreat.model)

m.w4.stereotype.spec <- boot.SE(data = wave4.new,
                                boots = 1000,
                                x = "lehDiff",
                                y = "stereotype",
                                pretreat.formula = pretreat.model.spec,
                                posttreat.formula = posttreat.model)

m.w5.stereotype.spec <- boot.SE(data = wave5.new,
                                boots = 1000,
                                x = "lehDiff",
                                y = "stereotype",
                                pretreat.formula = pretreat.model.spec,
                                posttreat.formula = posttreat.model)

# Export regression tables for SI (Table S2d.5)
stargazer(m.w4.allDVs.spec$model,
          m.w4.stereotype.spec$model,
          m.w4.malePref.spec$model,
          m.w5.allDVs.spec$model,
          m.w5.stereotype.spec$model,
          m.w5.malePref.spec$model,
          style = "ajps",
          se = list(m.w4.allDVs.spec$se,
                    m.w4.stereotype.spec$se,
                    m.w4.malePref.spec$se,
                    m.w5.allDVs.spec$se,
                    m.w5.stereotype.spec$se,
                    m.w5.malePref.spec$se),
          title = "Wave 4 and 5 Results (Alternative Measure of Settlement Type)",
          model.numbers = FALSE,
          digits = 3,
          omit.stat = c("f", "ser", "aic", "wald", 
                        "rsq", "adj.rsq"),
          star.char = c("*", "**"),
          star.cutoffs = c(0.05, 0.01))

##### ROBUSTNESS CHECKS 7: ALTERNATIVE MEASURE OF CHURCH EXPOSURE ##### 
###### Set new pretreatment model (w Church exposure from Schulz et al 2019) ###### 
posttreat.model.church  <- "sex+age+unemployed+incomeDecile+married+religiosity+educ+ChurchExp"

###### Sequential-g models for DVs (w alternative Church exp measure) ###### 
m.w4.allDVs.church <- boot.SE(data = wave4.Church,
                              boots = 1000,
                              x = "lehDiff",
                              y = "allDVs",
                              pretreat.formula = pretreat.model,
                              posttreat.formula = posttreat.model.church)

m.w5.allDVs.church <- boot.SE(data = wave5.Church,
                              boots = 1000,
                              x = "lehDiff",
                              y = "allDVs",
                              pretreat.formula = pretreat.model,
                              posttreat.formula = posttreat.model.church)

m.w4.malePref.church <- boot.SE(data = wave4.Church,
                                boots = 1000,
                                x = "lehDiff",
                                y = "malePreference",
                                pretreat.formula = pretreat.model,
                                posttreat.formula = posttreat.model.church)

m.w5.malePref.church <- boot.SE(data = wave5.Church,
                                boots = 1000,
                                x = "lehDiff",
                                y = "malePreference",
                                pretreat.formula = pretreat.model,
                                posttreat.formula = posttreat.model.church)

m.w4.stereotype.church <- boot.SE(data = wave4.Church,
                                  boots = 1000,
                                  x = "lehDiff",
                                  y = "stereotype",
                                  pretreat.formula = pretreat.model,
                                  posttreat.formula = posttreat.model.church)

m.w5.stereotype.church <- boot.SE(data = wave5.Church,
                                  boots = 1000,
                                  x = "lehDiff",
                                  y = "stereotype",
                                  pretreat.formula = pretreat.model,
                                  posttreat.formula = posttreat.model.church)

# Export regression tables for SI (Table S2d.6)
stargazer(m.w4.allDVs.church$model,
          m.w4.stereotype.church$model,
          m.w4.malePref.church$model,
          m.w5.allDVs.church$model,
          m.w5.stereotype.church$model,
          m.w5.malePref.church$model,
          style = "ajps",
          se = list(m.w4.allDVs.church$se,
                    m.w4.stereotype.church$se,
                    m.w4.malePref.church$se,
                    m.w5.allDVs.church$se,
                    m.w5.stereotype.church$se,
                    m.w5.malePref.church$se),
          title = "Wave 4 and 5 Results (Alternative Measure of Church Exposure)",
          model.numbers = FALSE,
          digits = 3,
          omit.stat = c("f", "ser", "aic", "wald", 
                        "rsq", "adj.rsq"),
          star.char = c("*", "**"),
          star.cutoffs = c(0.05, 0.01))

##### ROBUTNESS CHECKS 8: PLACEBO OUTCOMES ##### 
###### EU attitudes ###### 
m.w4.EUenlarge <- boot.SE(data = wave4.new,
                          boots = 1000,
                          x = "lehDiff",
                          y = "EUenlarge",
                          pretreat.formula = pretreat.model,
                          posttreat.formula = posttreat.model)

m.w5.EUenlarge <- boot.SE(data = wave5.new,
                          boots = 1000,
                          x = "lehDiff",
                          y = "EUenlarge",
                          pretreat.formula = pretreat.model,
                          posttreat.formula = posttreat.model)

m.w4.EUconfidence <- boot.SE(data = wave4.new,
                             boots = 1000,
                             x = "lehDiff",
                             y = "EUconfidence",
                             pretreat.formula = pretreat.model,
                             posttreat.formula = posttreat.model)

m.w5.EUconfidence <- boot.SE(data = wave5.new,
                             boots = 1000,
                             x = "lehDiff",
                             y = "EUconfidence",
                             pretreat.formula = pretreat.model,
                             posttreat.formula = posttreat.model)

# Export regression tables for SI (Table S2d.7)
stargazer(m.w4.EUenlarge$model,
          m.w5.EUenlarge$model,
          m.w4.EUconfidence$model,
          m.w5.EUconfidence$model,
          style = "ajps",
          se = list(m.w4.EUenlarge$se,
                    m.w5.EUenlarge$se,
                    m.w4.EUconfidence$se,
                    m.w5.EUconfidence$se),
          title = "Wave 4 and 5 Placebo Results (Attitudes Toward EU)",
          model.numbers = FALSE,
          digits = 3,
          omit.stat = c("f", "ser", "aic", "wald", 
                        "rsq", "adj.rsq"),
          star.char = c("*", "**"),
          star.cutoffs = c(0.05, 0.01))

###### Economic ideology ###### 
m.w4.BizOwn <- boot.SE(data = wave4.new,
                       boots = 1000,
                       x = "lehDiff",
                       y = "BizOwn",
                       pretreat.formula = pretreat.model,
                       posttreat.formula = posttreat.model)

m.w5.BizOwn <- boot.SE(data = wave5.new,
                       boots = 1000,
                       x = "lehDiff",
                       y = "BizOwn",
                       pretreat.formula = pretreat.model,
                       posttreat.formula = posttreat.model)

m.w4.ResProvide <- boot.SE(data = wave4.new,
                           boots = 1000,
                           x = "lehDiff",
                           y = "ResProvide",
                           pretreat.formula = pretreat.model,
                           posttreat.formula = posttreat.model)

m.w5.ResProvide <- boot.SE(data = wave5.new,
                           boots = 1000,
                           x = "lehDiff",
                           y = "ResProvide",
                           pretreat.formula = pretreat.model,
                           posttreat.formula = posttreat.model)

# Export regression tables for SI (Table S2d.8)
stargazer(m.w4.BizOwn$model,
          m.w5.BizOwn$model,
          m.w4.ResProvide$model,
          m.w5.ResProvide$model,
          style = "ajps",
          se = list(m.w4.BizOwn$se,
                    m.w5.BizOwn$se,
                    m.w4.ResProvide$se,
                    m.w5.ResProvide$se),
          title = "Wave 4 and 5 Placebo Results (Economic Ideology)",
          model.numbers = FALSE,
          digits = 3,
          omit.stat = c("f", "ser", "aic", "wald", 
                        "rsq", "adj.rsq"),
          star.char = c("*", "**"),
          star.cutoffs = c(0.05, 0.01))

##### ADDITIONAL TESTS 1: ONLY NATIVES SAMPLE ##### 
m.w4.allDVs_noImm <- boot.SE(data = filter(wave4.new,
                                           native == 1),
                             boots = 1000,
                             x = "lehDiff",
                             y = "allDVs",
                             pretreat.formula = pretreat.model,
                             posttreat.formula = posttreat.model)

m.w5.allDVs_noImm <- boot.SE(data = filter(wave5.new,
                                           native == 1),
                             boots = 1000,
                             x = "lehDiff",
                             y = "allDVs",
                             pretreat.formula = pretreat.model,
                             posttreat.formula = posttreat.model)

m.w4.malePref_noImm <- boot.SE(data = filter(wave4.new,
                                             native == 1),
                               boots = 1000,
                               x = "lehDiff",
                               y = "malePreference",
                               pretreat.formula = pretreat.model,
                               posttreat.formula = posttreat.model)

m.w5.malePref_noImm <- boot.SE(data = filter(wave5.new,
                                             native == 1),
                               boots = 1000,
                               x = "lehDiff",
                               y = "malePreference",
                               pretreat.formula = pretreat.model,
                               posttreat.formula = posttreat.model)

m.w4.stereotype_noImm <- boot.SE(data = filter(wave4.new,
                                               native == 1),
                                 boots = 1000,
                                 x = "lehDiff",
                                 y = "stereotype",
                                 pretreat.formula = pretreat.model,
                                 posttreat.formula = posttreat.model)

m.w5.stereotype_noImm <- boot.SE(data = filter(wave5.new,
                                               native == 1),
                                 boots = 1000,
                                 x = "lehDiff",
                                 y = "stereotype",
                                 pretreat.formula = pretreat.model,
                                 posttreat.formula = posttreat.model)

# exclude immigrants from empirical analyses (Table S3a.1)
stargazer(m.w4.allDVs_noImm$model,
          m.w4.stereotype_noImm$model,
          m.w4.malePref_noImm$model,
          m.w5.allDVs_noImm$model,
          m.w5.stereotype_noImm$model,
          m.w5.malePref_noImm$model,
          style = "ajps",
          se = list(m.w4.allDVs_noImm$se,
                    m.w4.stereotype_noImm$se,
                    m.w4.malePref_noImm$se,
                    m.w5.allDVs_noImm$se,
                    m.w5.stereotype_noImm$se,
                    m.w5.malePref_noImm$se),
          title = "Wave 4 and 5 Results (Exclude Immigrants)",
          model.numbers = FALSE,
          digits = 3,
          omit.stat = c("f", "ser", "aic", "wald", 
                        "rsq", "adj.rsq"),
          star.char = c("*", "**"),
          star.cutoffs = c(0.05, 0.01)) # full results for SI

##### ADDITIONAL TESTS 2: ONLY IMMIGRANTS SAMPLE ##### 
m.w4.allDVs_Imm <- boot.SE(data = filter(wave4.new,
                                         native == 0),
                           boots = 1000,
                           x = "lehDiff",
                           y = "allDVs",
                           pretreat.formula = pretreat.model,
                           posttreat.formula = posttreat.model)

m.w5.allDVs_Imm <- boot.SE(data = filter(wave5.new,
                                         native == 0),
                           boots = 1000,
                           x = "lehDiff",
                           y = "allDVs",
                           pretreat.formula = pretreat.model,
                           posttreat.formula = posttreat.model)

m.w4.malePref_Imm <- boot.SE(data = filter(wave4.new,
                                           native == 0),
                             boots = 1000,
                             x = "lehDiff",
                             y = "malePreference",
                             pretreat.formula = pretreat.model,
                             posttreat.formula = posttreat.model)

m.w5.malePref_Imm <- boot.SE(data = filter(wave5.new,
                                           native == 0),
                             boots = 1000,
                             x = "lehDiff",
                             y = "malePreference",
                             pretreat.formula = pretreat.model,
                             posttreat.formula = posttreat.model)

m.w4.stereotype_Imm <- boot.SE(data = filter(wave4.new,
                                             native == 0),
                               boots = 1000,
                               x = "lehDiff",
                               y = "stereotype",
                               pretreat.formula = pretreat.model,
                               posttreat.formula = posttreat.model)

m.w5.stereotype_Imm <- boot.SE(data = filter(wave5.new,
                                             native == 0),
                               boots = 1000,
                               x = "lehDiff",
                               y = "stereotype",
                               pretreat.formula = pretreat.model,
                               posttreat.formula = posttreat.model)

# only immigrants from empirical analyses (Table S3a.2)
stargazer(m.w4.allDVs_Imm$model,
          m.w4.stereotype_Imm$model,
          m.w4.malePref_Imm$model,
          m.w5.allDVs_Imm$model,
          m.w5.stereotype_Imm$model,
          m.w5.malePref_Imm$model,
          style = "ajps",
          se = list(m.w4.allDVs_Imm$se,
                    m.w4.stereotype_Imm$se,
                    m.w4.malePref_Imm$se,
                    m.w5.allDVs_Imm$se,
                    m.w5.stereotype_Imm$se,
                    m.w5.malePref_Imm$se),
          title = "Wave 4 and 5 Results (Only Immigrants)",
          model.numbers = FALSE,
          digits = 3,
          omit.stat = c("f", "ser", "aic", "wald", 
                        "rsq", "adj.rsq"),
          star.char = c("*", "**"),
          star.cutoffs = c(0.05, 0.01)) # full results for SI

# plot coefficient estimates for main paper (Figure 3)
coef_out_Imm.df <- data.frame(waves = rep(c("wave4",
                                            "wave5"),
                                          each = 3),
                              outcome = rep(c("allDVs",
                                              "stereotype",
                                              "malePreference"),
                                            2),
                              coef = NA,
                              se = NA,
                              lb = NA,
                              ub = NA,
                              imm = "imm")

coef_out_Imm.df$coef[1] <- coef(m.w4.allDVs_Imm$model)["lehDiff"]
coef_out_Imm.df$coef[2] <- coef(m.w4.stereotype_Imm$model)["lehDiff"]
coef_out_Imm.df$coef[3] <- coef(m.w4.malePref_Imm$model)["lehDiff"]
coef_out_Imm.df$coef[4] <- coef(m.w5.allDVs_Imm$model)["lehDiff"]
coef_out_Imm.df$coef[5] <- coef(m.w5.stereotype_Imm$model)["lehDiff"]
coef_out_Imm.df$coef[6] <- coef(m.w5.malePref_Imm$model)["lehDiff"]

coef_out_Imm.df$se[1] <- m.w4.allDVs_Imm$se["lehDiff"]
coef_out_Imm.df$se[2] <- m.w4.stereotype_Imm$se["lehDiff"]
coef_out_Imm.df$se[3] <- m.w4.malePref_Imm$se["lehDiff"]
coef_out_Imm.df$se[4] <- m.w5.allDVs_Imm$se["lehDiff"]
coef_out_Imm.df$se[5] <- m.w5.stereotype_Imm$se["lehDiff"]
coef_out_Imm.df$se[6] <- m.w5.malePref_Imm$se["lehDiff"]

coef_out_Imm.df$lb <- coef_out_Imm.df$coef - 1.96*coef_out_Imm.df$se
coef_out_Imm.df$ub <- coef_out_Imm.df$coef + 1.96*coef_out_Imm.df$se

coef_out_Imm.df$outcome <- factor(coef_out_Imm.df$outcome,
                                  levels = rev(c("allDVs",
                                                 "stereotype",
                                                 "malePreference")))

coef_out_noImm.df <- data.frame(waves = rep(c("wave4",
                                              "wave5"),
                                            each = 3),
                                outcome = rep(c("allDVs",
                                                "stereotype",
                                                "malePreference"),
                                              2),
                                coef = NA,
                                se = NA,
                                lb = NA,
                                ub = NA,
                                imm = "non-imm")

coef_out_noImm.df$coef[1] <- coef(m.w4.allDVs_noImm$model)["lehDiff"]
coef_out_noImm.df$coef[2] <- coef(m.w4.stereotype_noImm$model)["lehDiff"]
coef_out_noImm.df$coef[3] <- coef(m.w4.malePref_noImm$model)["lehDiff"]
coef_out_noImm.df$coef[4] <- coef(m.w5.allDVs_noImm$model)["lehDiff"]
coef_out_noImm.df$coef[5] <- coef(m.w5.stereotype_noImm$model)["lehDiff"]
coef_out_noImm.df$coef[6] <- coef(m.w5.malePref_noImm$model)["lehDiff"]

coef_out_noImm.df$se[1] <- m.w4.allDVs_noImm$se["lehDiff"]
coef_out_noImm.df$se[2] <- m.w4.stereotype_noImm$se["lehDiff"]
coef_out_noImm.df$se[3] <- m.w4.malePref_noImm$se["lehDiff"]
coef_out_noImm.df$se[4] <- m.w5.allDVs_noImm$se["lehDiff"]
coef_out_noImm.df$se[5] <- m.w5.stereotype_noImm$se["lehDiff"]
coef_out_noImm.df$se[6] <- m.w5.malePref_noImm$se["lehDiff"]

coef_out_noImm.df$lb <- coef_out_noImm.df$coef - 1.96*coef_out_noImm.df$se
coef_out_noImm.df$ub <- coef_out_noImm.df$coef + 1.96*coef_out_noImm.df$se

coef_out_noImm.df$outcome <- factor(coef_out_noImm.df$outcome,
                                    levels = rev(c("allDVs",
                                                   "stereotype",
                                                   "malePreference")))

coef_out_Imm_noImm <- rbind.data.frame(coef_out_noImm.df,
                                       coef_out_Imm.df)

coef_out_Imm.p <- ggplot(coef_out_Imm_noImm, aes(x = coef, y = outcome)) +
  geom_point(size = 0.8) +
  geom_errorbarh(aes(xmin = lb, xmax = ub),
                 height = 0) +
  facet_wrap(waves~imm, scale="fixed",
             labeller = as_labeller(c("wave4" = "Wave 4",
                                      "wave5" = "Wave 5",
                                      "imm" = "Immigrants",
                                      "non-imm" = "Natives"))) +
  geom_vline(xintercept = 0,
             linetype = "dashed") +
  scale_y_discrete(name = "",
                   labels = c("Pro-Female\nPreference",
                              "Gender Roles",
                              "Modern\nPro-Female Bias")) +
  scale_x_continuous(name = "Effect of Historical Pro-Female Bias",
                     limits = c(-0.85,0.7)) +
  coord_cartesian(ylim = c(1, 3), clip="off") +
  theme_bw() +
  theme(text=element_text(size=6),
        plot.margin = unit(c(1,1,1,1), "lines"))

# Export coef plot for comparison between immigrants & non-immigrants (Figure 3)
pdf("plots/coef_out_Imm.pdf", width = 3, height = 3)
coef_out_Imm.p
dev.off()

##### ADDITIONAL TESTS 3: BLACK DEATH SEVERITY ##### 
boot.SE.int <- function(data, # specify data to be analyzed on
                        boots, # number of iterations 
                        x, # specify predictor variable (character variable),
                        mod, # specify moderator (character variable)
                        y, # specify outcome variable (character variable)
                        pretreat.formula, # specify pretreat formula (RHS only, in character class)
                        posttreat.formula) # specify posttreat formula (RHS only, in character class)
{ 
  posttreat <- unlist(str_split(posttreat.formula, "\\+")) # create vector of posttreat variable names
  fs_main <- lm(as.formula(paste(y, 
                                 "~",
                                 x, "*", mod,
                                 "+ as.factor(country) +",
                                 pretreat.formula,
                                 "+",
                                 posttreat.formula)),
                data = data,
                weights=1/data$nObs, 
                na.action=na.exclude)
  demed.y_main <- data[,y] - as.matrix(data[,posttreat]) %*% coef(fs_main)[posttreat] # get demediated outcome
  data$demed.y_main <- as.numeric(unlist(demed.y_main)) # add demed y to databoot
  ss_main <- lm(as.formula(paste("demed.y_main",
                                 "~",
                                 x, "*", mod,
                                 "+",
                                 pretreat.formula)), # run second stage regression
                data = data,
                weights=1/data$nObs, 
                na.action=na.exclude)
  gboot <- matrix(NA, nrow=boots, ncol= length(names(coef(ss_main)))) # empty matrix to store bootstrapped values
  cov <- rep(NA, boots)
  set.seed(12345)
  for (b in 1:boots){
    databoot <- data[sample(1:nrow(data), replace = TRUE),]
    fs <- lm(as.formula(paste(y, 
                              "~",
                              x, "*", mod,
                              "+ as.factor(country) +",
                              pretreat.formula,
                              "+",
                              posttreat.formula)), # run first stage regression
             data = databoot,
             weights=1/databoot$nObs, 
             na.action=na.exclude)
    demed.y <- databoot[,y] - as.matrix(databoot[,posttreat]) %*% coef(fs)[posttreat] # get demediated outcome
    databoot$demed.y <- as.numeric(unlist(demed.y)) # add demed y to databoot
    ss <- lm(as.formula(paste("demed.y",
                              "~",
                              x, "*", mod,
                              "+",
                              pretreat.formula)), # run second stage regression
             data = databoot,
             weights=1/databoot$nObs, 
             na.action=na.exclude)
    gboot[b,] <- coef(ss) # store bootstrapped coefs
    #cov[b] <- cov(ss)[x,mod]
  }
  colnames(gboot) <- names(coef(ss))
  coef.mean <- apply(gboot,2,mean)
  return(list(model = ss_main,
              se = apply(gboot,2,sd),
              cov = crossprod((gboot[,x] - coef.mean[x]),
                              (gboot[,paste(x,mod,sep=":")] - coef.mean[paste(x,mod,sep=":")]))/(nrow(gboot)-1))) # return bootstrapped SEs
}

m.w4.allDVs_Black <- boot.SE.int(data = wave4.new,
                                 boots = 1000,
                                 x = "lehDiff",
                                 mod = "plagueSeverity",
                                 y = "allDVs",
                                 pretreat.formula = pretreat.model,
                                 posttreat.formula = posttreat.model)

m.w5.allDVs_Black <- boot.SE.int(data = wave5.new,
                                 boots = 1000,
                                 x = "lehDiff",
                                 mod = "plagueSeverity",
                                 y = "allDVs",
                                 pretreat.formula = pretreat.model,
                                 posttreat.formula = posttreat.model)

# Export regression tables for SI (Table S3b.1)
stargazer(m.w4.allDVs_Black$model,
          m.w5.allDVs_Black$model,
          style = "ajps",
          se = list(m.w4.allDVs_Black$se,
                    m.w5.allDVs_Black$se),
          title = "Black Death Full Results",
          model.numbers = FALSE,
          digits = 3,
          omit.stat = c("f", "ser", "aic", "wald", 
                        "rsq", "adj.rsq"),
          star.char = c("*", "**"),
          star.cutoffs = c(0.05, 0.01))

# Export coef plot for conditional marginal effects (Figure 4)
sim.data <- data.frame(sim = rep(seq(0.01,1,0.01),2),
                       waves = rep(c("wave4", "wave5"), each = 100))
sim.data$margins[c(1:100)] <- coef(m.w4.allDVs_Black$model)["lehDiff"] +
  coef(m.w4.allDVs_Black$model)["lehDiff:plagueSeverity"] * sim.data$sim[c(1:100)]
sim.data$margins[c(101:200)] <- coef(m.w5.allDVs_Black$model)["lehDiff"] +
  coef(m.w5.allDVs_Black$model)["lehDiff:plagueSeverity"] * sim.data$sim[c(101:200)]
sim.data$se <- NA
sim.data$se[c(1:100)] <- sqrt(m.w4.allDVs_Black$se["lehDiff"]^2 + 
                                m.w4.allDVs_Black$se["lehDiff:plagueSeverity"]^2 * (sim.data$sim[c(1:100)])^2 +
                                2 * sim.data$sim[c(1:100)] * as.numeric(m.w4.allDVs_Black$cov))
sim.data$se[c(101:200)] <- sqrt(m.w5.allDVs_Black$se["lehDiff"]^2 + 
                                  m.w5.allDVs_Black$se["lehDiff:plagueSeverity"]^2 * (sim.data$sim[c(101:200)])^2 +
                                  2 * sim.data$sim[c(1:100)] * as.numeric(m.w5.allDVs_Black$cov))
sim.data$lb <- sim.data$margins - 1.96*sim.data$se
sim.data$ub <- sim.data$margins + 1.96*sim.data$se

g.allDVs.Black.out.p <- ggplot(sim.data, aes(x = sim, 
                                             y = margins)) +
  geom_line() +
  geom_ribbon(aes(ymin = lb, ymax = ub),
              alpha = 0.3) +
  geom_hline(yintercept = 0,
             linetype = "dashed") +
  facet_wrap("waves", scales = "free_x",
             labeller = as_labeller(c("wave4" = "Wave 4",
                                      "wave5" = "Wave 5"))) +
  scale_x_continuous(name = "Black Death Severity") +
  scale_y_continuous(name = "Marginal Effect of\nHistorical Pro-Female Bias",
                     limits = c(-0.55,0.4)) +
  coord_cartesian(ylim = c(-0.55,0.4), clip="off") +
  theme_bw() +
  theme(text=element_text(size=6),
        plot.margin = unit(c(1,1,1,1), "lines"))

# Export coef plot for black death as pdf (Figure 4)
pdf("plots/g.allDVs.Black.out.pdf", width = 3, height = 2)
g.allDVs.Black.out.p
dev.off()

##### Only analyze sites that are closer (<= 160 km) ##### 
m.w4.allDVs_Black.short <- boot.SE.int(data = filter(wave4.new,
                                                     Dist_site_pollen <= 160),
                                       boots = 1000,
                                       x = "lehDiff",
                                       mod = "plagueSeverity",
                                       y = "allDVs",
                                       pretreat.formula = pretreat.model,
                                       posttreat.formula = posttreat.model)

m.w5.allDVs_Black.short <- boot.SE.int(data = filter(wave5.new,
                                                     Dist_site_pollen <= 160),
                                       boots = 1000,
                                       x = "lehDiff",
                                       mod = "plagueSeverity",
                                       y = "allDVs",
                                       pretreat.formula = pretreat.model,
                                       posttreat.formula = posttreat.model)

# Export coef plot for conditional marginal effects (SI)
sim.data2 <- data.frame(sim = rep(seq(0.01,1,0.01),2),
                        waves = rep(c("wave4", "wave5"), each = 100))
sim.data2$margins[c(1:100)] <- coef(m.w4.allDVs_Black.short$model)["lehDiff"] +
  coef(m.w4.allDVs_Black.short$model)["lehDiff:plagueSeverity"] * sim.data2$sim[c(1:100)]
sim.data2$margins[c(101:200)] <- coef(m.w5.allDVs_Black.short$model)["lehDiff"] +
  coef(m.w5.allDVs_Black.short$model)["lehDiff:plagueSeverity"] * sim.data2$sim[c(101:200)]
sim.data2$se <- NA
sim.data2$se[c(1:100)] <- sqrt(m.w4.allDVs_Black.short$se["lehDiff"]^2 + 
                                 m.w4.allDVs_Black.short$se["lehDiff:plagueSeverity"]^2 * (sim.data2$sim[c(1:100)])^2 +
                                2 * sim.data2$sim[c(1:100)] * as.numeric(m.w4.allDVs_Black.short$cov))
sim.data2$se[c(101:200)] <- sqrt(m.w5.allDVs_Black.short$se["lehDiff"]^2 + 
                                   m.w5.allDVs_Black.short$se["lehDiff:plagueSeverity"]^2 * (sim.data2$sim[c(101:200)])^2 +
                                  2 * sim.data2$sim[c(1:100)] * as.numeric(m.w5.allDVs_Black.short$cov))
sim.data2$lb <- sim.data2$margins - 1.96*sim.data2$se
sim.data2$ub <- sim.data2$margins + 1.96*sim.data2$se

g.allDVs.Black.out.short.p <- ggplot(sim.data2, aes(x = sim, 
                                                    y = margins)) +
  geom_line() +
  geom_ribbon(aes(ymin = lb, ymax = ub),
              alpha = 0.3) +
  geom_hline(yintercept = 0,
             linetype = "dashed") +
  facet_wrap("waves", scales = "free_x",
             labeller = as_labeller(c("wave4" = "Wave 4",
                                      "wave5" = "Wave 5"))) +
  scale_x_continuous(name = "Black Death Severity",
                     limits = c(0,0.5)) +
  scale_y_continuous(name = "Marginal Effect of\nHistorical Pro-Female Bias",
                     limits = c(-1,1)) +
  coord_cartesian(ylim = c(-1,1), clip="off") +
  theme_bw() +
  theme(text=element_text(size=6),
        plot.margin = unit(c(1,1,1,1), "lines"))

# Export coef plot for black death (closer sites) as pdf (Figure S3b.1)
pdf("plots/g.allDVs.Black.out.short.pdf", width = 3, height = 2)
g.allDVs.Black.out.short.p
dev.off()

##### Control for crop suitability ##### 
m.w4.allDVs_Black.crop <- boot.SE.int(data = wave4.new,
                                       boots = 1000,
                                       x = "lehDiff",
                                       mod = "plagueSeverity",
                                       y = "allDVs",
                                       pretreat.formula = pretreat.model.crop,
                                       posttreat.formula = posttreat.model)

m.w5.allDVs_Black.crop <- boot.SE.int(data = wave5.new,
                                       boots = 1000,
                                       x = "lehDiff",
                                       mod = "plagueSeverity",
                                       y = "allDVs",
                                       pretreat.formula = pretreat.model.crop,
                                       posttreat.formula = posttreat.model)

# Export coef plot for conditional marginal effects (SI)
sim.data3 <- data.frame(sim = rep(seq(0.01,1,0.01),2),
                        waves = rep(c("wave4", "wave5"), each = 100))
sim.data3$margins[c(1:100)] <- coef(m.w4.allDVs_Black.crop$model)["lehDiff"] +
  coef(m.w4.allDVs_Black.crop$model)["lehDiff:plagueSeverity"] * sim.data3$sim[c(1:100)]
sim.data3$margins[c(101:200)] <- coef(m.w5.allDVs_Black.crop$model)["lehDiff"] +
  coef(m.w5.allDVs_Black.crop$model)["lehDiff:plagueSeverity"] * sim.data3$sim[c(101:200)]
sim.data3$se <- NA
sim.data3$se[c(1:100)] <- sqrt(m.w4.allDVs_Black.crop$se["lehDiff"]^2 + 
                                 m.w4.allDVs_Black.crop$se["lehDiff:plagueSeverity"]^2 * (sim.data2$sim[c(1:100)])^2 +
                                 2 * sim.data3$sim[c(1:100)] * as.numeric(m.w4.allDVs_Black.crop$cov))
sim.data3$se[c(101:200)] <- sqrt(m.w5.allDVs_Black.crop$se["lehDiff"]^2 + 
                                   m.w5.allDVs_Black.crop$se["lehDiff:plagueSeverity"]^2 * (sim.data2$sim[c(101:200)])^2 +
                                   2 * sim.data3$sim[c(1:100)] * as.numeric(m.w5.allDVs_Black.crop$cov))
sim.data3$lb <- sim.data3$margins - 1.96*sim.data3$se
sim.data3$ub <- sim.data3$margins + 1.96*sim.data3$se

g.allDVs.Black.out.crop.p <- ggplot(sim.data3, aes(x = sim, 
                                                    y = margins)) +
  geom_line() +
  geom_ribbon(aes(ymin = lb, ymax = ub),
              alpha = 0.3) +
  geom_hline(yintercept = 0,
             linetype = "dashed") +
  facet_wrap("waves", scales = "free_x",
             labeller = as_labeller(c("wave4" = "Wave 4",
                                      "wave5" = "Wave 5"))) +
  scale_x_continuous(name = "Black Death Severity",
                     limits = c(0,0.38)) +
  scale_y_continuous(name = "Marginal Effect of\nHistorical Pro-Female Bias",
                     limits = c(-1,1)) +
  coord_cartesian(ylim = c(-1,1), clip="off") +
  theme_bw() +
  theme(text=element_text(size=6),
        plot.margin = unit(c(1,1,1,1), "lines"))

# Export coef plot for black death (with crop suitability) as pdf (Figure S3b.2)
pdf("plots/g.allDVs.Black.out.crop.pdf", width = 3, height = 2)
g.allDVs.Black.out.crop.p
dev.off()

##### ADDITIONAL TESTS 4: GENDER NORM PERSISTENCE IN THE AMERICAS ##### 
load("DataAnalysis_Americas.RData")

# recode education variable
placebo$educ <- ifelse(placebo$V248 == 1, 0,
                       ifelse(placebo$V248 == 2, 1,
                              ifelse(placebo$V248 == 3, 2,
                                     ifelse(placebo$V248 %in% c(4,6), 3,
                                            ifelse(placebo$V248 %in% c(5,7), 4,
                                                   ifelse(placebo$V248 == 8, 5,
                                                          ifelse(placebo$V248 == 9, 6, NA)))))))


table(placebo$settType)
placebo$Rural <- ifelse(placebo$settType == "Rural", 1,0)
placebo$Urban_large <- ifelse(placebo$settType == "Urban, large", 1, 0)
placebo$Urban_small <- ifelse(placebo$settType == "Urban, small", 1, 0)

# summary statistics (WVS, Table S3c.3)
sum_stats_wvs <- placebo %>%
  dplyr::select(allDVs,stereotype,malePreference,
                lehDiff,lowElevation,topography.x,
                postChristian,medieval,postMedieval,
                Rural,Urban_small,Urban_large,
                sex,age,unemployed,incomeDecile,
                married,religiosity,educ) %>%
  psych::describe() %>%
  data.frame() # summary statistics
sum_stats_wvs <- sum_stats_wvs[,c(2:5, 8:9,11)]
rownames(sum_stats_wvs) <- c("Modern Pro-Female Bias",
                             "Gender Roles",
                             "Female Priority",
                             "Historical Pro-Female Bias",
                             "Low Elevation",
                             "Topography",
                             "Post Christian",
                             "Medieval",
                             "Post-medieval",
                             "Rural",
                             "Urban (small)",
                             "Urban (large)",
                             "Sex",
                             "Age",
                             "Unemployed",
                             "Income",
                             "Married",
                             "Religiosity",
                             "Education")
xtable(sum_stats_wvs)

# run sequential-g regression models (WVS)
pretreat.model3 <- "lowElevation+as.factor(topography.x)+postChristian+medieval+postMedieval+as.factor(settType)"
posttreat.model3 <- "sex+age+unemployed+incomeDecile+married+religiosity+educ"

m.wvs.allDVs <- boot.SE(data = placebo,
                        boots = 1000,
                        x = "lehDiff",
                        y = "allDVs",
                        pretreat.formula = pretreat.model3,
                        posttreat.formula = posttreat.model3)

m.wvs.malePref <- boot.SE(data = placebo,
                          boots = 1000,
                          x = "lehDiff",
                          y = "malePreference",
                          pretreat.formula = pretreat.model3,
                          posttreat.formula = posttreat.model3)

m.wvs.stereotype <- boot.SE(data = placebo,
                            boots = 1000,
                            x = "lehDiff",
                            y = "stereotype",
                            pretreat.formula = pretreat.model3,
                            posttreat.formula = posttreat.model3)

# Export regression tables for SI (Table S3c.4)
stargazer(m.wvs.allDVs$model,
          m.wvs.stereotype$model,
          m.wvs.malePref$model,
          style = "ajps",
          se = list(m.wvs.allDVs$se,
                    m.wvs.stereotype$se,
                    m.wvs.malePref$se),
          title = "Western Bones Full Results",
          model.numbers = FALSE,
          digits = 3,
          omit.stat = c("f", "ser", "aic", "wald", 
                        "rsq", "adj.rsq"),
          star.char = c("*", "**"),
          star.cutoffs = c(0.05, 0.01))




